Advanced dimension reduction methods (including Principal Curves and surfaces, Independent Component Analysis, Non-negative Matrix Factorization)

Introduction

library(dplyr)

Attaching package: ‘dplyr’

The following object is masked from ‘package:Biobase’:

    combine

The following objects are masked from ‘package:BiocGenerics’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(ggplot2)
library(plotly)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(analogue) # Principal Curves wrapper

Loading required package: vegan
Loading required package: permute
Loading required package: lattice
This is vegan 2.5-7
analogue version 0.17-5

Attaching package: ‘analogue’

The following object is masked from ‘package:NMF’:

    compare

library(fastICA)
library(NMF)

data("iris")

Principal Curves

Fitting principal curves

# https://fromthebottomoftheheap.net/2014/01/09/pcurve-part-2/

x <- iris[1:4]
y <- iris[5]

results_iris_pc <- prcurve(x, plotit = TRUE, maxit = 100)

results_iris_pc

    Principal Curve Fitting

Call: prcurve(X = x, maxit = 100, plotit = TRUE)

Algorithm converged after 11 iterations

           SumSq Proportion
Total     681.37      1.000
Explained 650.45      0.955
Residual   30.92      0.045

Fitted curve uses 17.551 degrees of freedom.

Results in 2D

results_iris_pc_s <- data.frame(results_iris_pc$s[results_iris_pc$tag, ])
results_iris_pc_s['y'] <- y

ggplot() + 
  geom_line(data = results_iris_pc_s, mapping = aes(x = Petal.Length, y = Petal.Width)) +
  geom_point(data = iris, mapping = aes(x = Petal.Length, y = Petal.Width, color = Species))

Results in 3D

fig <- iris %>% 
  mutate(color = case_when(
    Species == 'setosa' ~ 'red',
    Species == 'virginica' ~ 'green',
    Species == 'versicolor' ~ 'blue'
  )) %>% 
  plot_ly(
    x = ~Sepal.Length,
    y = ~Sepal.Width,
    z = ~Petal.Width,
    type = 'scatter3d',
    mode = 'markers',
    marker = list(size = 1, color = ~color)
  ) 


results_iris_pc_smooth_curve <- results_iris_pc_s %>% mutate(color = 'black')
fig <- fig %>% 
  add_trace(
    data = results_iris_pc_smooth_curve,
    x = ~Sepal.Length,
    y = ~Sepal.Width,
    z = ~Petal.Width,
    type = 'scatter3d',
    mode = 'lines'
  )

fig

`arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.A marker object has been specified, but markers is not in the mode
Adding markers to the mode...
A marker object has been specified, but markers is not in the mode
Adding markers to the mode...

Results in Principal Space

plot(results_iris_pc)

# Fitted line in PC dimension
iris_s_pc <- predict(results_iris_pc$ordination, results_iris_pc$s, type = "wa", scaling = 0)[, 1:2] %>% as.data.frame()

iris_pc <- predict(results_iris_pc$ordination, x, type = "wa", scaling = 0)[, 1:2] %>% as.data.frame()
iris_pc['y'] <- y

fig <- ggplot() +
  geom_line(data = iris_s_pc, mapping = aes(x = PC1, y = PC2)) +
  geom_point(data = iris_pc,  mapping = aes(x = PC1, y = PC2, color=y))

fig

Principal Surface

There is no publicly available package.

ICA

Simple example

# http://rstudio-pubs-static.s3.amazonaws.com/93614_be30df613b2a4707b3e5a1a62f631d19.html
# https://rdrr.io/cran/fastICA/man/fastICA.html

# Source matrix
S <- cbind(sin((1:1000)/20), rep((((1:200)-100)/100), 5))
# Mixing matrix
A <- matrix(c(0.29, 0.655, -0.543, 0.557), 2, 2)
# plot graphs
par(mfcol = c(1, 2))
plot(1:1000, S[,1], type = "l",xlab = "S1", ylab = "")
plot(1:1000, S[,2], type = "l", xlab = "S2", ylab = "")

# Mixed two signals
X <- S %*% A

par(mfcol = c(1, 2))
plot(1:1000, X[,1], type = "l",xlab = "X1", ylab = "")
plot(1:1000, X[,2], type = "l", xlab = "X2", ylab = "")

# ICA for extracting independent sources from mixed signals
a <- fastICA(X, 2, alg.typ = "parallel", fun = "logcosh", alpha = 1,
             method = "R", row.norm = FALSE, maxit = 200,
             tol = 0.0001, verbose = FALSE)

par(mfcol = c(1, 2))
plot(1:1000, a$S[,1], type = "l", xlab = "S'1", ylab = "")
plot(1:1000, a$S[,2], type = "l", xlab = "S'2", ylab = "")

par(mfcol = c(2, 3))
plot(1:1000, S[,1 ], type = "l", main = "Original Signals", 
     xlab = "", ylab = "")

Iris example

x <- iris[1:4]
y <- iris[5]

results_iris_ica <- fastICA(x, n.comp = 2)

A list containing the following components:

• Pre-whitening matrix that projects data onto the first n.comp principal components.

results_iris_ica$K

           [,1]       [,2]
[1,] -0.1763375 -1.3373258
[2,]  0.0412425 -1.4871770
[3,] -0.4180098  0.3531217
[4,] -0.1748261  0.1537381

• Estimated un-mixing matrix

results_iris_ica$W

           [,1]      [,2]
[1,] -0.9940472 0.1089499
[2,]  0.1089499 0.9940472

• Estimated mixing matrix

results_iris_ica$A

           [,1]       [,2]       [,3]        [,4]
[1,]  0.7010963 -0.2112469  1.7544864  0.73394561
[2,] -0.4011386 -0.3374820 -0.1066651 -0.04316124

• Estimated source matrix

S <- data.frame(results_iris_ica$S)
S['class'] <- y

ggplot(S, aes(x = X1, y = X2, color = class)) + 
  geom_point()

NMF

Iris example

x <- iris[1:4]
y <- iris[5]

results_iris_nmf <- nmf(x, rank = 2)
results_iris_nmf

<Object of class: NMFfit>
 # Model:
  <Object of class:NMFstd>
  features: 150 
  basis/rank: 2 
  samples: 4 
 # Details:
  algorithm:  brunet 
  seed:  random 
  RNG: 10403L, 355L, ..., 1268043011L [4975d7ab2202c0cdbd26424e1365bcc4]
  distance metric:  'KL' 
  residuals:  3.086644 
  Iterations: 430 
  Timing:
     user  system elapsed 
     0.03    0.00    0.03 

results_iris_nmf_H <- results_iris_nmf@fit@H
results_iris_nmf_H

     Sepal.Length Sepal.Width Petal.Length Petal.Width
[1,]    1.0401571  0.72681648    0.2563836  0.02780619
[2,]    0.4111733  0.08199243    0.5652751  0.21979063

results_iris_nmf_W <- data.frame(results_iris_nmf@fit@W)
results_iris_nmf_W

results_iris_nmf_W['y'] <- y

ggplot() +
  geom_point(data = results_iris_nmf_W, aes(x = X1, y = X2, color = y))